{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np \n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preliminaries\n",
    "\n",
    "We first load up the merged data from SHARE and HRS. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_pickle('../data_sources/hrs-share_wide_select.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the function for a weighted mean. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wmean(x,var,wvar):\n",
    "    xx = x.loc[~x[var].isna(),:]\n",
    "    names = {var: (xx[var] * xx[wvar]).sum()/xx[wvar].sum()}\n",
    "    return pd.Series(names, index=[var])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We prepare a simple table with stats, the 2-year transition rates and the fraction in good health in wave 1 and 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnames = ['DE','SE','NL','SP','IT','FR','DK','US']\n",
    "table = pd.DataFrame(index=cnames,columns=['gg2','gb2','g_w1','g_w2'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute weighted means by country for these vars. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "for c in table.columns:\n",
    "\ttable[c] = df.groupby('cname').apply(wmean,var=c,wvar='wgid_w1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute one-year transition rates\n",
    "\n",
    "We use the eigenvalues and vectors for the two year transition matrix to compute the one-year transition matrix. We also compute the steady state. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "for c in table.index:\n",
    "\tp2 = np.zeros((2,2))\n",
    "\tp2[0,0] = table.loc[c,'gg2']\n",
    "\tp2[0,1] = 1-table.loc[c,'gg2']\n",
    "\tp2[1,0] = table.loc[c,'gb2']\n",
    "\tp2[1,1] = 1- p2[1,0]\n",
    "\teigvalues,eigvectors = np.linalg.eig(p2)\n",
    "\teigvalues = np.sqrt(eigvalues)\n",
    "\tp1 = eigvectors @ np.diag(eigvalues) @ np.linalg.inv(eigvectors)\n",
    "\tps = np.linalg.matrix_power(p1,1000)\n",
    "\ttable.loc[c,'gg'] = p1[0,0]\n",
    "\ttable.loc[c,'gb'] = p1[1,0]\n",
    "\ttable.loc[c,'g_s'] = ps[0,0] \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the following stats for one year transition rates and implied steady-state fraction in good health (g_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gg2</th>\n",
       "      <th>gb2</th>\n",
       "      <th>g_w1</th>\n",
       "      <th>g_w2</th>\n",
       "      <th>gg</th>\n",
       "      <th>gb</th>\n",
       "      <th>g_s</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>DE</th>\n",
       "      <td>0.957591</td>\n",
       "      <td>0.533106</td>\n",
       "      <td>0.933130</td>\n",
       "      <td>0.928282</td>\n",
       "      <td>0.974321</td>\n",
       "      <td>0.322796</td>\n",
       "      <td>0.926311</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SE</th>\n",
       "      <td>0.954206</td>\n",
       "      <td>0.500441</td>\n",
       "      <td>0.945582</td>\n",
       "      <td>0.931213</td>\n",
       "      <td>0.972638</td>\n",
       "      <td>0.299017</td>\n",
       "      <td>0.916165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NL</th>\n",
       "      <td>0.960598</td>\n",
       "      <td>0.712649</td>\n",
       "      <td>0.948468</td>\n",
       "      <td>0.950634</td>\n",
       "      <td>0.973696</td>\n",
       "      <td>0.475751</td>\n",
       "      <td>0.947608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SP</th>\n",
       "      <td>0.943659</td>\n",
       "      <td>0.568025</td>\n",
       "      <td>0.930858</td>\n",
       "      <td>0.915159</td>\n",
       "      <td>0.965069</td>\n",
       "      <td>0.352178</td>\n",
       "      <td>0.909763</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT</th>\n",
       "      <td>0.953888</td>\n",
       "      <td>0.527967</td>\n",
       "      <td>0.924979</td>\n",
       "      <td>0.923286</td>\n",
       "      <td>0.972098</td>\n",
       "      <td>0.319472</td>\n",
       "      <td>0.919676</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FR</th>\n",
       "      <td>0.956852</td>\n",
       "      <td>0.521163</td>\n",
       "      <td>0.927794</td>\n",
       "      <td>0.925796</td>\n",
       "      <td>0.974008</td>\n",
       "      <td>0.313941</td>\n",
       "      <td>0.923538</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DK</th>\n",
       "      <td>0.967526</td>\n",
       "      <td>0.569275</td>\n",
       "      <td>0.930578</td>\n",
       "      <td>0.942490</td>\n",
       "      <td>0.980090</td>\n",
       "      <td>0.349019</td>\n",
       "      <td>0.946034</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>US</th>\n",
       "      <td>0.932617</td>\n",
       "      <td>0.360515</td>\n",
       "      <td>0.891752</td>\n",
       "      <td>0.870705</td>\n",
       "      <td>0.961635</td>\n",
       "      <td>0.205261</td>\n",
       "      <td>0.842526</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         gg2       gb2      g_w1      g_w2        gg        gb       g_s\n",
       "DE  0.957591  0.533106  0.933130  0.928282  0.974321  0.322796  0.926311\n",
       "SE  0.954206  0.500441  0.945582  0.931213  0.972638  0.299017  0.916165\n",
       "NL  0.960598  0.712649  0.948468  0.950634  0.973696  0.475751  0.947608\n",
       "SP  0.943659  0.568025  0.930858  0.915159  0.965069  0.352178  0.909763\n",
       "IT  0.953888  0.527967  0.924979  0.923286  0.972098  0.319472  0.919676\n",
       "FR  0.956852  0.521163  0.927794  0.925796  0.974008  0.313941  0.923538\n",
       "DK  0.967526  0.569275  0.930578  0.942490  0.980090  0.349019  0.946034\n",
       "US  0.932617  0.360515  0.891752  0.870705  0.961635  0.205261  0.842526"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We check the assumption of stationarity by looking at whether the fraction in good health in steady-state is close to the fraction in good health in the data. The ranking is quite good. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeM0lEQVR4nO3df5gdVZ3n8feHmJgOSWgwLW4SIMGJgQgC2psRdRlWRhIRJDKsgyuzY0CZDL/iqhnFGUZdZwfcqA88woIREBGUGTVEZNGIjIKOP0iHTggJRGMASYeRIITwI0oSvvtHVcPl2n3v6b5dfet2f17Pc5++darq1vekO/3tOufUOYoIzMzM6tmr2QGYmVlrcMIwM7MkThhmZpbECcPMzJI4YZiZWZKXNTuAoTRlypSYMWNGs8MwM2sZq1evfiwiOlKOHVEJY8aMGXR1dTU7DDOzliHpodRj3SRlZmZJnDDMzCyJE4aZmSVxwjAzsyROGGZmlmREjZIyMxtNVnT3sHTlRrZu38nU9jaWzJvNgqOmFXY9Jwwzsxa0oruHC5avY+euPQD0bN/JBcvXARSWNNwkZWbWgpau3PhCsui1c9celq7cWNg1fYdhZoM23E0i9qKt23cOqHwo+A7DzAalt0mkZ/tOghebRFZ09zQ7tFFhanvbgMqHghOGmQ1KM5pE7EVL5s2mbeyYl5S1jR3DknmzC7umm6TMbFCa0SRiL+pt+vMoKTMrvantbfT0kRyKbBKxl1pw1LRh7TNyk5SZDUozmkSsuXyHYWaD0owmEWsuJwwzG7ThbhKx5nKTlJmZJXHCMDOzJE4YZmaWpNCEIWm+pI2SNkn6WB/795V0k6R7JN0l6bCq/WMkdUu6pcg4zcysvsIShqQxwOXA24E5wHskzak67OPAmoh4HfA/gEur9i8G7isqRjMzS1fkHcZcYFNEbI6I54AbgZOrjpkD3A4QEfcDMyTtDyBpOvAO4KoCYzQzs0RFJoxpwMMV21vyskprgVMAJM0FDgKm5/suAf4OeL7AGM3MLFGRCUN9lEXV9sXAvpLWAOcB3cBuSScCj0bE6roXkc6S1CWpa9u2bY3GbGZm/Sjywb0twAEV29OBrZUHRMQOYCGAJAEP5K/TgHdKOgEYD0yWdH1EnF59kYhYBiwD6OzsrE5IZmY2RIq8w1gFzJI0U9I4siRwc+UBktrzfQDvB+6MiB0RcUFETI+IGfl5/9ZXsjAzs+FT2B1GROyWdC6wEhgDXBMR6yUtyvdfCRwKXCdpD7ABOLOoeMzMrDGKGDmtOJ2dndHV1dXsMMzMWoak1RHRmXKsn/Q2M7MkThhmZpbECcPMzJI4YZiZWZK6CUPSKZJ+JelJSTskPSVpx3AEZ2Zm5ZEyrPb/ACdFhCcBNDMbxVKapH7rZGFmZv3eYUg6JX/bJelfgBXAH3r3R8TyYkMzM7MyqdUkdVLF+2eB4yu2A3DCMDMbRfpNGBHROyngmyPi3yv3SXpz0YGZmVm5pPRhfCGxzMzMRrBafRhHA28COiR9qGLXZLLJBM3MbBSp1YcxDpiYHzOponwHcGqRQZmZWfnU6sO4A7hD0rUR8dAwxmRmZiVUq0nqO+RLqmaL4b1URLyzuLDMzKxsajVJfXbYojAzs9Kr1yRlZmYGJMwlJWkWcBEwBxjfWx4RBxcYl5mNYCu6e1i6ciNbt+9kansbS+bNZsFR05odltWR8hzGl4ErgN3AfwWuA75aZFBmNnKt6O7hguXr6Nm+kwB6tu/kguXrWNHd0+zQrI6UhNEWEbeTrf/9UER8EnhrsWGZ2Ui1dOVGdu7a85Kynbv2sHTlxiZFZKlSpjf/vaS9gF9JOhfoAV5ZbFhmNlJt3b5zQOVWHil3GB8EJgDnA28ATgf+usCYzGwEm9reNqByK4+6CSMiVkXE08ATEbEwIv4iIn4+DLGZ2Qi0ZN5s2sa+dHahtrFjWDJvdpMislQpo6SOBq4mmybkQElHAH8TEWcXHZyZFW+4Ryz1frZHSbWelD6MS4B5wM0AEbFW0jFFBmVmw6N3xFJvJ3TviCWg8KThBNF6UvowiIiHq4r29HmgmbUUj1iygUi5w3hY0puAkDSOrPPba3ybjQAesWQDkXKHsQg4B5gGbAGOzLfNrMV5xJINRMooqcci4r0RsX9EvDIiTo+I3w1HcGZWLI9YsoFIGSXVAXwAmFF5fESckXDufOBSshX6roqIi6v27wtcA7wa+D1wRkTcK+kAsilIXgU8DyyLiEsT62RmiTxiyQYipQ/j28CPgR8wgM5uSWOAy4G3kTVlrZJ0c0RsqDjs48CaiHiXpEPy448jm7fqwxFxt6RJwGpJt1Wda2a5RobGesSSpUpJGBMi4qOD+Oy5wKaI2Awg6UbgZKDyl/4csplwiYj7Jc2QtH9EPAI8kpc/Jek+sj4UJwyzKs0aGmujT0qn9y2SThjEZ08DKofjbsnLKq0FTgGQNBc4CJheeYCkGcBRwC/6uoiksyR1Seratm3bIMI0a20eGmvDpd+EIekpSTuAxWRJY6ekHRXl9fzxuq75kq8VLgb2lbQGOA/oJmuO6o1hIvAt4IMR0ec1I2JZRHRGRGdHR0dCWGYji4fG2nCpteLepAY/ewtwQMX2dGBr1TV2AAsBlC0c/kD+QtJYsmRxQ0QsbzAWsxFransbPX0kBw+NtaGW9KT3IK0CZkmamT/wdxr59CK9JLXn+wDeD9wZETvy5HE1cF9EfL7AGM1anofG2nBJ6fQelIjYna+fsZJsWO01EbFe0qJ8/5XAocB1kvaQdWifmZ/+ZuCvgHV5cxXAxyPi1qLiNWtVHhrb2lppuVpFVHcrtK7Ozs7o6upqdhhmZkmqR7hBdnd40SmHD1vSkLQ6IjpTjq3V6b1frdfQhWtmNjq12gi3Wk1Sq8lGNQk4EHgif98O/AaYWXRwZmYjWauNcOv3DiMiZkbEwWR9ECdFxJSIeAVwIuBRS2ZmDWq1yR9TRkn958rO5oj4LvBnxYVkZjY6tNoIt5RRUo9J+gfgerImqtMBz1ZrZtagVhvhlpIw3gN8Argp374zLzMzswa10uSPdRNGRDwOLJY0GXg+Ip4uPiwzMyubun0Ykg6X1A2sA9ZLWi3psOJDMzOzMknp9P4i8KGIOCgiDgI+DCwrNiwzMyublISxd0T8sHcjIn4E7F1YRGZmVkopnd6bJV0IfDXfPp18RlkzM2it+ZBs8FLuMM4AOsge1rspf7+wyKDMrHX0zofUs30nwYsr/q3o7ml2aDbEUkZJPQGc71FSZtaXWvMh+S5jZPEoKTNrSKvNh2SD51FSZtaQVpsPyQbPo6TMrCGtNh+SDZ5HSZlZQ1ptPiQbvJSEcQbwKbJRUiKbS8qjpMzsBa00H5INXvIoqWGIxczMSqxuwpD0HbJpzSs9CXQBX4yI3xcRmJmZlUtKp/dm4GngS/lrB/Bb4DX5tpmZjQIpfRhHRcQxFdvfkXRnRBwjaX1RgZmZWbmk3GF0SDqwdyN/PyXffK6QqMzMrHRS7jA+DPxE0q/JRknNBM6WtDfwlSKDMzOz8kgZJXWrpFnAIWQJ4/6Kju5LCozNzMxKJOUOg4j4A7C24FjMzKzEUvowzMzMnDDMzCxNUpOUpGnAQZXHR8SdRQVlZmblk/Kk92eAvwQ2AL2rpATZnFL1zp0PXAqMAa6KiIur9u8LXAO8Gvg9cEZE3JtyrpmZDa+UO4wFwOy84zuZpDHA5cDbgC3AKkk3R8SGisM+DqyJiHdJOiQ//rjEc83MbBilTg0ydhCfPRfYFBGbI+I54Ebg5Kpj5gC3A0TE/cAMSfsnnmtmZsMo5Q7jWWCNpNuBF+4yIqLeDLbTgIcrtrcAf1p1zFrgFLIHA+eS9ZNMTzwXAElnAWcBHHjggX0dYmZmQyAlYdycvwZKfZRVz3p7MXCppDVka4Z3A7sTz80KI5aRLxnb2dnZ5zFmZta4lCe9vyJpHNnstAAbI2JXwmdvAQ6o2J4ObK367B3kizFJEtlKfg8AE+qda2Zmw6tuH4akY4FfkXVC/1/gl5KOqXVObhUwS9LMPOGcRtWdiqT2fB/A+4E78yRS91wzMxteKU1SnwOOj4iNAJJeA3wdeEOtkyJit6RzgZVkQ2OviYj1khbl+68EDgWuk7SHbNjumbXOHUwFzcxsaKQkjLG9yQIgIn4pKWnUVETcCtxaVXZlxfufAbNSzzUzs+ZJSRhdkq4GvppvvxdYXVxIZmZWRikJ42+Bc4DzyUYv3UnWl2FmZqNIyiipPwCfz19mZjZK9ZswJP1rRLxb0jr6eAYiIl5XaGRmZlYqte4wFudfTxyOQMzMrNz6fQ4jIh7J354dEQ9VvoCzhyc8MzMri5TJB9/WR9nbhzoQMzMrt1p9GH9Ldifxakn3VOyaBPy06MDMzKxcavVhfA34LnAR8LGK8qci4vFCozIzs9Kp1YfxZEQ8SLbq3eMV/Re7JPU51biZmY1cKX0YVwBPV2w/k5eZmdkokpIwFBEvPIcREc+T9oS4mZmNIElLtEo6X9LY/LWYbNlWMzMbRVISxiLgTUAPLy6VelaRQZmZWfmkzCX1KNkCRmZmNorVTRiSxpMtbPRaYHxveUScUWBcZmZWMilNUl8FXgXMA+4gW1/7qSKDMjOz8klJGH8SERcCz0TEV4B3AIcXG5aZmZVNSsLYlX/dLukwYB9gRmERmZlZKaU8T7FM0r7APwA3AxOBCwuNyszMSqdmwpC0F7AjIp4gW5r14GGJyszMSqdmk1T+VPe5wxSLmZmVWEofxm2SPiLpAEn79b4Kj8zMzEolpQ+j93mLcyrKAjdPmZmNKrUWUPpvEfEN4LiI8NxRZmajXK0mqQvyr98cjkDMzKzcajVJ/U7SD4GZkm6u3hkR7ywuLDMzK5taCeMdwOvJpgb53PCEY2ZmZdVvwoiI54CfS3pTRGwbzIdLmk+2xOsY4KqIuLhq/z7A9cCBeSyfjYgv5/v+J/B+sg72dcDCiPj9YOIwM7PG1R1W20CyGANcDrwdmAO8R9KcqsPOATZExBHAscDnJI2TNA04H+iMiMPIEo6nWDcza6KU5zAGay6wKSI253crNwInVx0TwCRJIpty5HFgd77vZUCbpJcBE4CtBcZqZmZ11E0YDTykNw14uGJ7S15W6TLgULJksA5YHBHPR0QP8FngN8AjwJMR8f1+4jtLUpekrm3bBnUzZGZmCVLuMH4h6RuSTsjvBFL1dWxUbc8D1gBTgSOByyRNzic7PBmYme/bW9LpfV0kIpZFRGdEdHZ0dAwgPDMzG4iUhPEaYBnwV8AmSf8s6TUJ520BDqjYns4fNystBJZHZhPwAHAI8OfAAxGxLSJ2AcvJ1hU3M7MmSen0joi4LSLeQzZq6a+BuyTdIenoGqeuAmZJmilpHFmndfXzHL8BjgOQtD8wG9icl79R0oT8ruY44L4B1s3MzIZQyprerwBOJ7vD+C1wHtkv/iOBb5A1G/2RiNgt6VxgJdkop2siYr2kRfn+K4FPA9dKWkfWhPXRiHgMeEzSN4G7yTrBu8nucszMrEkUUd2tUHWA9Euyh/e+HBFbqvZ9NCI+U2B8A9LZ2RldXV3NDsPMrGVIWh0RnSnH1ltAaQxwS0R8uq/9ZUoWZmZWrHoLKO0BjhimWMzMrMRS1sNYk08++A3gmd7CiFheWFRmZlY6KQljP+B3wFsryoJsqKuZmY0SdRNGRCwcjkDMzKzcUobVjgfOBF4LjO8tj4gz+j3JzMxGnJQnvb8KvIpsGo87yJ7YfqrIoMzMrHxSEsafRMSFwDMR8RWyhZUOLzYsMzMrm5SEsSv/ul3SYcA+wIzCIjIzs1JKGSW1LJ899kKyKUEmAv9YaFRmZlY6KaOkrsrf3gEcXGw4ZmZWVimjpF4O/AVZM9QLx0fE/youLDMzK5uUJqlvA08Cq4E/FBuOmZmVVUrCmB4R8wuPxMzMSi1llNRPJXkYrZnZKJdyh/EW4H2SHiBrkhLZQnyvKzQyMzMrlZSE8fbCozAzs9JLWdP7IaAdOCl/tedlZmY2itRNGJIWAzcAr8xf10s6r+jAzMysXFKapM4E/jQingGQ9BngZ8AXigzMzMzKJWWUlIA9Fdt78jIzMxtFUu4wvgz8QtJN+fYC4OrCIjIzs1JKmUvq85J+RDa8VsDCiOguOjAzMyuXfhOGpMkRsUPSfsCD+at3334R8Xjx4ZmZWVnUusP4GnAi2RxSUVGufNsz15qZjSL9JoyIODH/OnP4wjEzs7JKeQ7j9pQyMzMb2Wr1YYwHJgBT8hX3eofSTgamDkNsZmZWIrX6MP4G+CBZcljNiwljB3B5sWGZmVnZ9NskFRGX5v0XH4mIgyNiZv46IiIuS/lwSfMlbZS0SdLH+ti/j6TvSForab2khRX72iV9U9L9ku6TdPSgamhmZkMi5TmML0g6DJgDjK8ov67WeZLGkN2JvA3YAqySdHNEbKg47BxgQ0ScJKkD2Cjphoh4DrgU+F5EnCppHFnzmJmZNUnKmt6fAI4lSxi3kk13/hOgZsIA5gKbImJz/jk3AicDlQkjgEmSBEwEHgd2S5oMHAO8DyBPIM+lVsrMzIZeylxSpwLHAf8REQuBI4CXJ5w3DXi4YntLXlbpMuBQYCuwDlgcEc+TPeOxDfiypG5JV0nau6+LSDpLUpekrm3btiWEZWZmg5GSMHbmv8R7//J/lLSH9vqaoDCqtucBa8g61o8ELsuv8TLg9cAVEXEU8AzwR30gABGxLCI6I6Kzo6MjISwzMxuMlITRJakd+BLZaKm7gbsSztsCHFCxPZ3sTqLSQmB5ZDYBDwCH5OduiYhf5Md9kyyBmJlZk6R0ep+dv71S0veAyRFxT8JnrwJmSZoJ9ACnAf+96pjfkDV3/VjS/sBsYHNEPCbpYUmzI2JjfswGzMysaVI6vW+PiOMAIuLB6rL+RMRuSecCK4ExwDURsV7Sonz/lcCngWslrSNrwvpoRDyWf8R5wA35CKnNZHcjZmbWJIU+6R0Rt5KNrKosu7Li/Vbg+H7OXQN0plzHzMyK5ye9zcwsSa3Zai8FLpV0XkR4/W4zs1EuZZTUf0iaBCDpHyQtl+QRS2Zmo0xKwrgwIp6S9Bay5ya+AlxRbFhmZlY2KQljT/71HWQP0n0bGFdcSGZmVkYpCaNH0heBdwO3Snp54nlmZjaCpPzifzfZsxTzI2I7sB+wpMigzMysfFKe9H4WWF6x/QjwSJFBmZlZ+bhpyczMkjhhmJlZEicMMzNL4oRhZmZJnDDMzCyJE4aZmSVxwjAzsyROGGZmlsQJw8zMkjhhmJlZEicMMzNL4oRhZmZJnDDMzCyJE4aZmSVxwjAzsyROGGZmlqTuAkoj3YruHpau3MjW7TuZ2t7GknmzWXDUtGaHZWZWOqM6Yazo7uGC5evYuWsPAD3bd3LB8nUAThpmZlVGdZPU0pUbX0gWvXbu2sPSlRubFJGZWXmN6oSxdfvOAZWbmY1mozphTG1vG1C5mdloVmjCkDRf0kZJmyR9rI/9+0j6jqS1ktZLWli1f4ykbkm3FBHfknmzaRs75iVlbWPHsGTe7CIuZ2bW0gpLGJLGAJcDbwfmAO+RNKfqsHOADRFxBHAs8DlJ4yr2LwbuKyrGBUdN46JTDmdaexsCprW3cdEph7vD28ysD0WOkpoLbIqIzQCSbgROBjZUHBPAJEkCJgKPA7vz46cD7wD+N/ChooJccNQ0JwgzswRFNklNAx6u2N6Sl1W6DDgU2AqsAxZHxPP5vkuAvwOepwZJZ0nqktS1bdu2oYjbzMz6UGTCUB9lUbU9D1gDTAWOBC6TNFnSicCjEbG63kUiYllEdEZEZ0dHR4Mhm5lZf4pMGFuAAyq2p5PdSVRaCCyPzCbgAeAQ4M3AOyU9CNwIvFXS9QXGamZmdRSZMFYBsyTNzDuyTwNurjrmN8BxAJL2B2YDmyPigoiYHhEz8vP+LSJOLzBWMzOro7BO74jYLelcYCUwBrgmItZLWpTvvxL4NHCtpHVkTVgfjYjHiorJzMwGTxHV3QqtS9I24KGq4inASExCrlfrGIl1AterldSq00ERkdQBPKISRl8kdUVEZ7PjGGquV+sYiXUC16uVDFWdRvXUIGZmls4Jw8zMkoyGhLGs2QEUxPVqHSOxTuB6tZIhqdOI78MwM7OhMRruMMzMbAg4YZiZWZKWThiNrLchabGke/PyDw5r4DUk1GlfSTdJukfSXZIOSz23mRqs1zWSHpV07/BGXd9g6yXpAEk/lHRf/jO4ePij71sDdRqfb/f+f/vU8Effv0Z+BvP9ha7PM1gN/t96UNI6SWskddW9WES05Ivs6fFfAwcD44C1wJyqYz4OfCZ/30E2ffo44DDgXmAC2dPuPwBmtUidlgKfyN8fAtyeem4r1ivfPgZ4PXBvs+syhN+v/wS8Pn8/CfhlGb5fDdZJwMT8/VjgF8Abm12nofgZzMs+BHwNuKXZ9RmqegEPAlNSr9fKdxgvrLcREc+RTVJ4ctUx/a23cSjw84h4NiJ2A3cA7xq+0PuVUqc5wO0AEXE/MCOfhyvl3GZppF5ExJ1k37uyGXS9IuKRiLg7L3+KbKGwMizM0kidIiKezo8Zm7/KMqqmoZ9Bvbg+z1XDF3KShuo1UK2cMBpZb+Ne4BhJr5A0ATiBl86s2ywpdVoLnAIgaS5wENlMwCnnNksj9SqzIamXpBnAUWR/kTdbQ3XKm23WAI8Ct0VEGeoEjX+vLiFhfZ4maLReAXxf0mpJZ9W7WCsnjEGvtxER9wGfAW4Dvkf2D7q7sEjTpdTpYmDf/D/leUA3Wewp5zZLI/Uqs4brJWki8C3ggxGxo6A4B6KhOkXEnog4kuwX0tzqfoAmGnS9NID1eZqg0Z/BN0fE68mW0j5H0jG1LlbkEq1FS11v4+LIGus2Sepdb+OuiLgauBpA0j/nn9dsdeuU/1JZCJA3tT2QvybUO7eJGqlXmTVUL0ljyZLFDRGxfDgCTjAk36uI2C7pR8B8sjv6ZmukXqeRrc9zAjAemCzp+ijHkgsNfb8iYmv+9VFJN5E1cd3Z79Wa3WnTQGfPy4DNwExe7Ox5bdUxVwCfzN/vD/SQd/AAr8y/HgjcD+zbInVqB8bl7z8AXJd6bivWq2L/DMrX6d3I90vAdcAlza7HENapA2jP37cBPwZObHadhupnMC8/lnJ1ejfy/dobmFTx/qfA/JrXa3aFG/zHOoFsdMmvgb/PyxYBi/L3U4Hvk/Vf3AucXnHuj4EN+T/wcc2uywDqdDTwK7Ikt5yKRNfXuWV5NVivrwOPALvI/qI6s9n1abRewFvImg7uIWs2XQOc0Oz6NFin15E1d9yT/3/7x2bXZah+Bis+41hKlDAa/H4dnP/+WwusT/md4alBzMwsSSt3epuZ2TBywjAzsyROGGZmlsQJw8zMkjhhmJlZEicMGzEkfVLSR+ocs0DSnOGKqZb+4pV0raRTh+gaD0qaIqld0tkV5ceWbdZVKz8nDBttFpBNxjbatANn1zvIrBYnDGtpkv4+XwvgB8DsivIPSFqVr83wLUkTJL0JeCewNJ///9V9HdfHNTok3SbpbklflPSQpCn5vg8pW1flXlWsq1KjvM94+3CMpJ9K2lx5tyFpSR7vPapYb0LSinwCufX9TCJ3MfDqvN5L87KJkr4p6X5JN+TTRpj1r9lPKfrl12BfwBvInuKfAEwGNgEfyfe9ouK4fwLOy99fC5xasa/P46qucxlwQf5+PtkT2lMqrr832fT568lmna1X/kfxVl3vWuAbZH/QzSGbvhrgeGAZ2bQiewG3AMfk+/bLv7aRPWX9inz7wTzWGVRMrUL2xPKTZHMP7QX8DHhLs7+nfpX71cqTD5r9F+CmiHgWQNLNFfsOk/RPZE0xE4GV/XxGynFvIV8vJSK+J+mJivKbIuKZ/PrL85jUT/leNeKttiKyqfg3VKxdcHz+6s63JwKzyCaLO19S75ouB+Tlv6vx+ZBNwrklj2UNWVL5SZ1zbBRzwrBW19/cNtcCCyJiraT3kf1FPdjj+muqGWg5pE85/4c+Pk/ARRHxxZdcTDoW+HPg6Ih4Np8ldvwAr7EH/z6wOtyHYa3sTuBdktokTQJOqtg3CXgkn0L8vRXlT+X76h1X6SfAuwEkHQ/sW3H9BXn/yN5kdyE/rlPeX7wpVgJn5GtoIGmapFcC+wBP5MniEOCNfZxbXW+zAfNfFNayIuJuSf9CNtPrQ2S/lHtdSLaC3UNk/Qa9vyxvBL4k6Xzg1BrHVfoU8HVJf0m2nO8jwFP59a8F7sqPuyoiuiEbGttPeX/xptT3+5IOBX6W908/DZxOtgjYIkn3ABuBn/dx7u8k/buke4HvAv9vINc2AzxbrVk9kl4O7ImI3ZKOBq6IbFU5s1HFdxhm9R0I/KukvYDnyBahMRt1fIdhZmZJ3OltZmZJnDDMzCyJE4aZmSVxwjAzsyROGGZmluT/A7wePhWPj+8bAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.scatter(table['g_w1'],table['g_s'])\n",
    "plt.xlabel('data good health')\n",
    "plt.ylabel('stationary fraction good health')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient\n",
    "\n",
    "For the static model in section 2, we need information on the gradient. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/loulou/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1667: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = value\n"
     ]
    }
   ],
   "source": [
    "for c in table.index:\n",
    "\tdf_c = df.loc[df['cname']==c,:]\n",
    "\tdf_c.loc[:,'qinc'] = pd.qcut(df_c.loc[:,'hitot_w2'],q=4)\n",
    "\tgrad = df_c.groupby('qinc').apply(wmean,var='g_w1',wvar='wgid_w1')\n",
    "\tgrad.index = [x for x in range(1,5)]\n",
    "\ttable.loc[c,'g_q2'] = grad.loc[2,'g_w1']/grad.loc[1,'g_w1']\n",
    "\ttable.loc[c,'g_q3'] = grad.loc[3,'g_w1']/grad.loc[1,'g_w1']\n",
    "\ttable.loc[c,'g_q4'] = grad.loc[4,'g_w1']/grad.loc[1,'g_w1']\n",
    "\tincs = df_c.groupby('qinc').apply(wmean,var='hitot_w1',wvar='wgid_w1')\n",
    "\tincs.index = [x for x in range(1,5)]\n",
    "\ttable.loc[c,'inc_q1'] = incs.loc[1,'hitot_w1']\n",
    "\ttable.loc[c,'inc_q2'] = incs.loc[2,'hitot_w1']\n",
    "\ttable.loc[c,'inc_q3'] = incs.loc[3,'hitot_w1']\n",
    "\ttable.loc[c,'inc_q4'] = incs.loc[4,'hitot_w1']\n",
    "\tmincs = incs['hitot_w1'].mean()\n",
    "\tfor q in range(1,5):\n",
    "\t\ttable.loc[c,'inc_q'+str(q)] *= 1/mincs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>g_q2</th>\n",
       "      <th>g_q3</th>\n",
       "      <th>g_q4</th>\n",
       "      <th>inc_q1</th>\n",
       "      <th>inc_q2</th>\n",
       "      <th>inc_q3</th>\n",
       "      <th>inc_q4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>DE</th>\n",
       "      <td>1.017634</td>\n",
       "      <td>1.061784</td>\n",
       "      <td>1.078631</td>\n",
       "      <td>0.561241</td>\n",
       "      <td>0.670246</td>\n",
       "      <td>0.963283</td>\n",
       "      <td>1.805230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SE</th>\n",
       "      <td>1.050808</td>\n",
       "      <td>1.061996</td>\n",
       "      <td>1.074089</td>\n",
       "      <td>0.542441</td>\n",
       "      <td>0.879170</td>\n",
       "      <td>1.114809</td>\n",
       "      <td>1.463580</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NL</th>\n",
       "      <td>1.030899</td>\n",
       "      <td>1.024509</td>\n",
       "      <td>1.023527</td>\n",
       "      <td>0.587189</td>\n",
       "      <td>0.820913</td>\n",
       "      <td>1.085506</td>\n",
       "      <td>1.506392</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SP</th>\n",
       "      <td>1.025151</td>\n",
       "      <td>0.997648</td>\n",
       "      <td>1.051461</td>\n",
       "      <td>0.585783</td>\n",
       "      <td>0.631921</td>\n",
       "      <td>1.183983</td>\n",
       "      <td>1.598313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT</th>\n",
       "      <td>1.025042</td>\n",
       "      <td>1.014742</td>\n",
       "      <td>1.038783</td>\n",
       "      <td>0.590575</td>\n",
       "      <td>0.859482</td>\n",
       "      <td>0.990178</td>\n",
       "      <td>1.559765</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FR</th>\n",
       "      <td>1.030915</td>\n",
       "      <td>1.060553</td>\n",
       "      <td>1.092331</td>\n",
       "      <td>0.509560</td>\n",
       "      <td>0.795265</td>\n",
       "      <td>1.124371</td>\n",
       "      <td>1.570804</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DK</th>\n",
       "      <td>1.087277</td>\n",
       "      <td>1.079270</td>\n",
       "      <td>1.101802</td>\n",
       "      <td>0.490115</td>\n",
       "      <td>0.799720</td>\n",
       "      <td>1.144787</td>\n",
       "      <td>1.565379</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>US</th>\n",
       "      <td>1.174827</td>\n",
       "      <td>1.223211</td>\n",
       "      <td>1.275940</td>\n",
       "      <td>0.338666</td>\n",
       "      <td>0.630960</td>\n",
       "      <td>0.955165</td>\n",
       "      <td>2.075210</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        g_q2      g_q3      g_q4    inc_q1    inc_q2    inc_q3    inc_q4\n",
       "DE  1.017634  1.061784  1.078631  0.561241  0.670246  0.963283  1.805230\n",
       "SE  1.050808  1.061996  1.074089  0.542441  0.879170  1.114809  1.463580\n",
       "NL  1.030899  1.024509  1.023527  0.587189  0.820913  1.085506  1.506392\n",
       "SP  1.025151  0.997648  1.051461  0.585783  0.631921  1.183983  1.598313\n",
       "IT  1.025042  1.014742  1.038783  0.590575  0.859482  0.990178  1.559765\n",
       "FR  1.030915  1.060553  1.092331  0.509560  0.795265  1.124371  1.570804\n",
       "DK  1.087277  1.079270  1.101802  0.490115  0.799720  1.144787  1.565379\n",
       "US  1.174827  1.223211  1.275940  0.338666  0.630960  0.955165  2.075210"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table[['g_q2','g_q3','g_q4','inc_q1','inc_q2','inc_q3','inc_q4']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing moments \n",
    "\n",
    "We compute the mean and standard deviation of moments using the bootstrap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def onesamp(df):\n",
    "\ttable = pd.DataFrame(index=cnames,columns=['gg2','gb2'])\n",
    "\t# sample with replacement\n",
    "\tdfi = df.groupby('cname', group_keys=False).apply(lambda x: x.sample(n=len(x),replace=True))\n",
    "\t# compute stats\n",
    "\tfor c in table.columns:\n",
    "\t\ttable[c] = dfi.groupby('cname').apply(wmean,var=c,wvar='wgid_w1')\n",
    "\tfor c in table.index:\n",
    "\t\tp2 = np.zeros((2,2))\n",
    "\t\tp2[0,0] = table.loc[c,'gg2']\n",
    "\t\tp2[0,1] = 1-p2[0,0]\n",
    "\t\tp2[1,0] = table.loc[c,'gb2']\n",
    "\t\tp2[1,1] = 1- p2[1,0]\n",
    "\t\teigvalues,eigvectors = np.linalg.eig(p2)\n",
    "\t\teigvalues = np.sqrt(eigvalues)\n",
    "\t\tp1 = eigvectors @ np.diag(eigvalues) @ np.linalg.inv(eigvectors)\n",
    "\t\ttable.loc[c,'gg'] = p1[0,0]\n",
    "\t\ttable.loc[c,'gb'] = p1[1,0]\n",
    "\t\tdf_c = dfi.loc[dfi['cname']==c,:]\n",
    "\t\tdf_c.loc[:,'qinc'] = pd.qcut(df_c.loc[:,'hitot_w2'],q=4)\n",
    "\t\tgrad = df_c.groupby('qinc').apply(wmean,var='g_w1',wvar='wgid_w1')\n",
    "\t\tgrad.index = [x for x in range(1,5)]\n",
    "\t\ttable.loc[c,'g_q2'] = grad.loc[2,'g_w1']/grad.loc[1,'g_w1']\n",
    "\t\ttable.loc[c,'g_q3'] = grad.loc[3,'g_w1']/grad.loc[1,'g_w1']\n",
    "\t\ttable.loc[c,'g_q4'] = grad.loc[4,'g_w1']/grad.loc[1,'g_w1']\n",
    "\ttable = table[['gg','gb','g_q2','g_q3','g_q4']]\n",
    "\treturn table\n",
    "\t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/loulou/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1667: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = value\n"
     ]
    }
   ],
   "source": [
    "moments = pd.DataFrame(columns=['rep','gg','gb','g_q2','g_q3','g_q4'])\n",
    "nreps = 1000\n",
    "np.random.seed(1234)\n",
    "for r in range(nreps):\n",
    "\tmom_r = onesamp(df)\n",
    "\tmom_r['rep'] = int(r)\n",
    "\tmoments = moments.append(mom_r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "moments = moments.reset_index()\n",
    "moments.set_index(['index','rep'],inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gg</th>\n",
       "      <th>gb</th>\n",
       "      <th>g_q2</th>\n",
       "      <th>g_q3</th>\n",
       "      <th>g_q4</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>index</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>DE</th>\n",
       "      <td>0.974090</td>\n",
       "      <td>0.327992</td>\n",
       "      <td>1.029659</td>\n",
       "      <td>1.071441</td>\n",
       "      <td>1.086319</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DK</th>\n",
       "      <td>0.979945</td>\n",
       "      <td>0.349289</td>\n",
       "      <td>1.082923</td>\n",
       "      <td>1.083623</td>\n",
       "      <td>1.102843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FR</th>\n",
       "      <td>0.974267</td>\n",
       "      <td>0.315457</td>\n",
       "      <td>1.031429</td>\n",
       "      <td>1.056177</td>\n",
       "      <td>1.090645</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT</th>\n",
       "      <td>0.971952</td>\n",
       "      <td>0.321932</td>\n",
       "      <td>1.022728</td>\n",
       "      <td>1.010093</td>\n",
       "      <td>1.039689</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NL</th>\n",
       "      <td>0.973459</td>\n",
       "      <td>0.476829</td>\n",
       "      <td>1.031051</td>\n",
       "      <td>1.025663</td>\n",
       "      <td>1.025642</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SE</th>\n",
       "      <td>0.972339</td>\n",
       "      <td>0.302608</td>\n",
       "      <td>1.052746</td>\n",
       "      <td>1.063812</td>\n",
       "      <td>1.076431</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SP</th>\n",
       "      <td>0.965110</td>\n",
       "      <td>0.356402</td>\n",
       "      <td>1.027981</td>\n",
       "      <td>1.000185</td>\n",
       "      <td>1.054194</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>US</th>\n",
       "      <td>0.961600</td>\n",
       "      <td>0.205977</td>\n",
       "      <td>1.177713</td>\n",
       "      <td>1.222356</td>\n",
       "      <td>1.276081</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             gg        gb      g_q2      g_q3      g_q4\n",
       "index                                                  \n",
       "DE     0.974090  0.327992  1.029659  1.071441  1.086319\n",
       "DK     0.979945  0.349289  1.082923  1.083623  1.102843\n",
       "FR     0.974267  0.315457  1.031429  1.056177  1.090645\n",
       "IT     0.971952  0.321932  1.022728  1.010093  1.039689\n",
       "NL     0.973459  0.476829  1.031051  1.025663  1.025642\n",
       "SE     0.972339  0.302608  1.052746  1.063812  1.076431\n",
       "SP     0.965110  0.356402  1.027981  1.000185  1.054194\n",
       "US     0.961600  0.205977  1.177713  1.222356  1.276081"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "moments.groupby('index').mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the standard deviations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gg</th>\n",
       "      <th>gb</th>\n",
       "      <th>g_q2</th>\n",
       "      <th>g_q3</th>\n",
       "      <th>g_q4</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>index</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>DE</th>\n",
       "      <td>0.003491</td>\n",
       "      <td>0.044133</td>\n",
       "      <td>0.030888</td>\n",
       "      <td>0.025817</td>\n",
       "      <td>0.025538</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DK</th>\n",
       "      <td>0.003571</td>\n",
       "      <td>0.048601</td>\n",
       "      <td>0.029608</td>\n",
       "      <td>0.028468</td>\n",
       "      <td>0.028135</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FR</th>\n",
       "      <td>0.003513</td>\n",
       "      <td>0.033703</td>\n",
       "      <td>0.023113</td>\n",
       "      <td>0.023630</td>\n",
       "      <td>0.020761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT</th>\n",
       "      <td>0.003938</td>\n",
       "      <td>0.044698</td>\n",
       "      <td>0.025370</td>\n",
       "      <td>0.025610</td>\n",
       "      <td>0.023636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NL</th>\n",
       "      <td>0.004240</td>\n",
       "      <td>0.068089</td>\n",
       "      <td>0.018193</td>\n",
       "      <td>0.018542</td>\n",
       "      <td>0.018212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SE</th>\n",
       "      <td>0.003707</td>\n",
       "      <td>0.042688</td>\n",
       "      <td>0.019814</td>\n",
       "      <td>0.021411</td>\n",
       "      <td>0.019188</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SP</th>\n",
       "      <td>0.005645</td>\n",
       "      <td>0.043020</td>\n",
       "      <td>0.029936</td>\n",
       "      <td>0.033832</td>\n",
       "      <td>0.028219</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>US</th>\n",
       "      <td>0.001506</td>\n",
       "      <td>0.009638</td>\n",
       "      <td>0.016617</td>\n",
       "      <td>0.016421</td>\n",
       "      <td>0.016628</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             gg        gb      g_q2      g_q3      g_q4\n",
       "index                                                  \n",
       "DE     0.003491  0.044133  0.030888  0.025817  0.025538\n",
       "DK     0.003571  0.048601  0.029608  0.028468  0.028135\n",
       "FR     0.003513  0.033703  0.023113  0.023630  0.020761\n",
       "IT     0.003938  0.044698  0.025370  0.025610  0.023636\n",
       "NL     0.004240  0.068089  0.018193  0.018542  0.018212\n",
       "SE     0.003707  0.042688  0.019814  0.021411  0.019188\n",
       "SP     0.005645  0.043020  0.029936  0.033832  0.028219\n",
       "US     0.001506  0.009638  0.016617  0.016421  0.016628"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "moments.groupby('index').std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We save these results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "moments.groupby('index').std().to_pickle('../estimation/moments/stds_health.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "moments.groupby('index').mean().to_pickle('../estimation/moments/means_health.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "cf2a50979671a58939829e6829efb726aa5da11149213b77bd50351f899d04fb"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
